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ABSTRACT 

The initial mass function (IMF), binary fraction and distributions of binary parameters (mass ratios, 
separations and eccentricities) are indispensable input for simulations of stellar populations. It is often 
claimed that these are poorly constrained significantly affecting evolutionary predictions. Recently, 
dedicated observing campaigns provided new constraints on the initial conditions for massive stars. 
Findings include a larger close binary fraction and a stronger preference for very tight systems. We 
investigate the impact on the predicted merger rates of neutron stars and black holes. Despite the 
changes with previous assumptions, we only find an increase of less than a factor 2 (insignificant 
compared with evolutionary uncertainties of typically a factor 10 — 100). We further show that the 
uncertainties in the new initial binary properties do not significantly affect (within a factor of 2) our 
predictions of double compact object merger rates. An exception is the uncertainty in IMF (variations 
by a factor of 6 up and down). No significant changes in the distributions of final component masses, 
mass ratios, chirp masses and delay times are found. We conclude that the predictions are, for practical 
purposes, robust against uncertainties in the initial conditions concerning binary parameters with 
exception of the IMF. This eliminates an important layer of the many uncertain assumptions affecting 
the predictions of merger detection rates with the gravitational wave detectors aLIGO/aVirgo. 

Subject headings: stars: black holes, neutron stars, x-ray binaries 


1. INTRODUCTION 

These are very exciting times for gravitational wave 
astrophysics. The direct detection of the gravitational 
wave signal of the merger of two compact objects, neu¬ 
tron stars (NS) or black holes (BH) is anticipated before 
the end of this decade. 

Gravitational waves are a natural consequ ence of the 
theory of General Relativity (Einstein 1918). They are 


perturbations of the spacetime metric propagating at the 
speed of light, which are generated, for example, during 
the inspiral of two compact objects. Indirect evidence 
for the existence of gravitational waves was pr ovided by 
the orbital decay of the Hulse-Taylor pulsar (Hulse & 


Taylor 1975 


Taylor & Weisberg 1989) a nd later with 


stronger constraints by the double pulsar (Burgay et al. 


2003 Lyne et al. 2004). Direct detection of the grav¬ 


itational wave signal of the inspiral of NS-NS, BH-NS 
or BH-BH binaries and the subsequent merger and ring 
down is expected to happen in the next few years, now 
the advanced ground-based gravitati onal wave detectors 
aLIGO and Virgo are coming online (Abbott et al.|2009 
Caron et al.||1997|). 


'The initial LICO/Virgo observations were concluded in 
2010 without detecti on, but they provid ed upper limits 
on the merger rates (Abadie et alj^012). The advanced 


version of the detectors will be approximately 10 times 
more sensitive than the initial versions, expanding the 
detection volume and thus the chance of dete ction by a 
factor of about a thousand (lAasi et al.||2013c I. 


The first sc ience run of adva nced LlCU is scheduled 
for late 2015 (Aasi et aL]|2013b). A detection during the 


first science run is not considered to be likely, but one or 
more detections are anticipated in the next few years as 
the sensitivity increases to a range of 200 Mpc for double 
neutron stars. This translates to ~ 0.2-200 exp ected 
detections d ouble neutron star mergers per year dAasif 


et al. 2013b). Even without detections, the new upper 
limits will become astrophysical interesting as they start 
to rul e out the models that predict the highest merger 
rates (Mandel &: O’Shaughn essy 2010 Belczynski et al. 
2012a Stevenson et al.||2015 ~ 


Obtaining reliable predictions of the merger rates of 
relativistic compact ob iects has been a ver y large chal¬ 
lenge, as reviewed by (Abadie et al. 2010). The rates 
quoted above are derived semi-empiricall y asing the ob¬ 


served binary pulsa rs in our Galaxy (Phinney 1991 
Kalogera et al. 20041. The large uncertainties result from 
the small size of the sample of observed binary pulsars 
and the pulsar luminosity distribution. The implicit as¬ 
sumption is made that the observed sample is represen¬ 
tative for the Galactic population. Further constraints 
come from short gamma-ray bursts, but the rate esti¬ 
mates depend on the uncertai n luminosity funct ion and 
opening angle of the jets (e.g. Fong et al. 2012). Type 
Ibc supernovae have served so far as ultimate upper lim¬ 
its (Kim et al.||2010). The merger rates for BH-NS and 
BH-BH systems are even more problematic as we lack di¬ 
rect observational evidence of their existence. However, 


some immediate progenitors for BH-BH (e.g. Bulik et aL | 
2011) and BH-NS systems (e.g. Grudzinska et al. 2015[ ) 


have been proposed. 

Even though theoretical predictions suffer from even 
larger uncertainties, they have been crucial to estimate 
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the merger rates involving black holes. They further pro¬ 
vide the expected distribution of properties for all merger 
types including for example the delay times, component 
and chirp masses. More importantly, the theoretical 
predictions are providing the tool for future compari¬ 
son with the detections, crucial to identify what grav¬ 
itational wave signals teach us about the astrophy sics of 
the progenitor systems (e.g. Stevenson et al.|2015|. Sev- 


in the past decade (e.g. [Lipunov et al. 

1997 

Bethe & 

Brown|19981 

De Donder &: Vanbeveren|1998 


loom et al. 


Tauris|2003[ Dewi fc Pols|2603[ Nutzman et al.|2604[|De 

Donder & Vanbeveren||2004t IPfahl et al. 1120051 IPostnov 


¥ 


Yungelson|2006l |Yungi 

al.||2008l 120101 IDomin 


elson et al.|2006||b’Shaughnessy 


et al. 1120081 |2010r bominik et al.||20121 2^13[ [Mennekens 

& Vanbeveren||2014| |Dominik et al.pOlS). 

hbr all predictions we can distinguish several sources 
of uncertainty: (I) the adopted initial conditions in the 
simulations, (II) the uncertainties in the physics of the 
stellar evolution and binary interaction, (III) the uncer¬ 
tainties associated with the normalization of the merger 
rates, which include the mapping to the appropriate star 
formation history of the detection volumes and (IV) mod¬ 
eling of the detectability of gravitational waves from the 
predicted merger events. 

Recent observations have provided new constraints on 
the initial conditions and primordial binary properties 


Ghmi et al.||2U12| [Sana et al. 

IlSana e( al.||ii01ij 

20141 IKobulnicky et al.|2014[|JVI 

oe & Di Stefano|2015a|b 


a review see 


raus 


2(J13p. i'he studies show that the initial conditions tor 


young massive stars differ substantially from the initial 
conditions that have typically been adopted to simulate 
compact object mergers coming from binary evolutionary 
channels. 

Among the most striking findings are (i) the large frac¬ 
tion massive stars that have a companion close enough 
to interact by exchanging mass before they die and (ii) 
the preference for very tight binaries with orbital periods 
of a few days, implying that a large fraction of massive 
stars will interact even before leaving the main sequence. 
Further findings include (hi) the confirmation of a flat 
distribution of mass ratios, ruling out a distribution that 
is strongly peaked to equal masses and (iv) a distribu¬ 
tion of eccentricities that favors circular systems in stark 
contrast with the typically adopted thermal eccentricity 
distribution which favors eccentric systems. 

The first two conclusions depend on adopted model 
of stellar evolution and in particular they require at 
least modest radial expansion for the majority of massive 
stars. The amount by which stars expand is consider¬ 
ably uncertain for high mass stars. Note that some very 
rapidly rotating m assive stars may actua lly decrease in 
size as they evolve ( Yoon fc Langer||2005 1, possibly even 


preventing interactio n with companion by mass transfer 
(de Mink et al.|2009 1 in very tight systems. For example, 


let us consider the fraction of massive binary systems in 
which mass transfer starts already during the main se¬ 
quence evolution of the primary star. This depends on 
the maximum radius that a star reaches on the main se¬ 
quence which depends on a number of rather uncertain 
processes such as convection and overshooting, (e.g. cali- 


bratio ns by |Pols et al.|I997] [Ribas et al.|2000||Brott et al.| 


by rotation (e.! 

r. Yoon & Langer|^2005l iBrott et al. 

|2011b 

Ekstrom et al.j 

2012 Szecsi et al. 12015), mass loss 

by stel- 


possible density inversions in the outer layers of massive 


stars which can result in infla ted envelopes(e.g. Yusof 
et al. 20I3[ [Kohler et al. 2015[ and Jiang et al. 21)15 


subm.). Varying the maximum radius of main sequence 
stars by 30% up and down and allowing for uncertainties 
in the initial distributions, we find a large variation for 
this fraction, 20% — 50%Q 

A further caveat to keep in mind is that the obser¬ 
vations are limited to regions nearby such as our own 
Galaxy and the Magellanic Clouds which only cover 
metallicities down to about one fifth of the solar value. 
Also the regime of the hig hest mass stars is not well 


probed. Sana et al. (2012) observations cover stars in 
mass range 15 — bUM©, but most of these are towards 
the lower end of this mass range. For low metallici¬ 
ties and higher star masses we have no direct constraints 
and the uncertainties in the evolutionary models start to 
play a larger role. Despite several attempts, no trends 
with metallicity or environment are found so far (e.g. 


Bastian et al. 


Moe & Di Stefano 2013). We ex¬ 


tend" Sana et al. (|2U12[) distributions in our study to the 


entire range that can produce double compact objects 
(Afzams = 5 — I 5 OM 0 ) and to a broad range of metallic¬ 
ities {Z = 0.002 — 0.02) as discussed in Section 2. 

In this study we focus on the binary formation chan- 


that may occur in dense star clusters, e.g. 

Ivanova et al. 

2008 

[Banerjee et al. 112010 [Aarseth 

20121 Clausen et al. 

2013 

jSamsing et al.|2014 Ramirez- 

luiz et 

al.|2015|). We 


of the new initial conditions? and (ii) How robust are 
the predictions against the allowed variations in the new 
initial conditions resulting from the observational uncer¬ 
tainties. 

For this purpose we perform a comparative p opulation 
synthesis stu dy where we use the recent work byjDominikj 


et al. (20121 as a reference. We simulate the evolution of 
massive biliary systems following the evolutionary chan¬ 
nels for the formation of double compact objects. In 
Sect. [ 2 ] we describe the new and old initial conditions. In 
Sect, [^we give a brief description of our computational 
method, the physical assumptions and computation of 
the merger rates. In Sect. we compare the old and new 
initial distributions for the entire simulated populations 
and for the progenitors of double compact object merg¬ 
ers. In Sect, [^we discuss the impact on the resulting 
merger rates. Finally, in Sect. we present our discus¬ 
sion and conclusions. 


2. INITIAL DISTRIBUTIONS 
2.1. New standard model (N) and its variations 

Recent dedicated observing campaigns have provided 
new constraints on the binary properties of young mas- 

^ This estimate was done with the StarTrack population synthe¬ 
sis code described in Section 3 considering systems with primary 
masses in the range 15 — ISOMq. The initial distributions of pe¬ 
riod and eccentricity were altered within limits listed in Section 2.1 
assuming a binary fraction of 100%. 
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sive stars. Here, we investigate th e impact of t he dis - 
tributions obtained by the work of Sana et al. (2012). 


This study is based on an intense spectroscopic moni¬ 
toring spanning a decade surveying all 0-type stars in 
six nearby (< 3 kpc) very young (about 2 Myr old) open 
star clusters and associations. Even though the sam¬ 
ple may seem of modest size it exceeds previous studies 
in completeness in terms of the fraction of systems for 
which orbital solutions have been obtained. It provided 
an average of 20 radial velocity measurements for all 71 
systems in the open clusters/associations that contain at 
least one 0-star. This dataset includes orbital solutions 
for several long period systems (between 100 and 1000 
days) which are very challenging as they require a long 
term observing campaign. This sample allowed a robust 
derivation of the underlying distribution of binary pa¬ 
rameters after correction for incompleteness and biases. 

The very young ages, relatively low densities and ve¬ 
locity dispersion of the stars, imply that the effects of 
stellar evolution and dynamical interactions are minimal. 
This makes this sample the most suitable to provide con¬ 
straints on the primordial binary properties and thus the 
initial conditions for our simulations. 

The primary stars in this sample have spectral types 
ranging from 09.7 to 03, which correspond approxi¬ 
mately to a mass range from 15 to 60 Mq. This is 
appropriate for the progenitor systems of double com¬ 
pact object mergers that involve at least one black hole. 
Lower mass binaries may dominate the formation of dou¬ 
ble neutron stars. 

The most suitabl e study for this mass range is provided 
by Dunstall et al. (2015) for the early B-type stars in 
th e 30 Uor regionand the C ygnus OB2 sample analyzed 


by Kobulnicky et al. 


(2014), which contains stars down 


to spectral type s of B2.5V, which a pproximately corre¬ 


sponds to 8 Mq. Wright et al. (20151 infer an age spread 
for this region of 1 -7 My r. S tatistically the hnding s by 
Kobulnicky et al. (2014) and Dunst all et al. (|20 15 ) ar e 
consistent with the distr ibutions by Sana et ak] ( 2012 ), 
although the results by Kobulnicky et al. ”(2014) and 
Dunstall et al. (2015) favor a hatter period distribution 
and a lower close biliary fraction, similar to the old ini¬ 
tial distributions that we use as reference (see Sec. 2.2). 
Whether this is a sign of a trend with decreasing primary 
mass, or whether the distribution of ages play a role is 
not clear. 

Sourcing from our conclusions, we hnd that the 
changes in the period distribution have a rather small 
impact on double compa ct merger rates. Thus we ap¬ 
ply the Sana et al.|(2012) distributions as our new initial 
conditions tor both U arid early B stars for consistency. 


Orbital periods — For the dist ribution of o r bital p eriods, 
p, we adopt the distribution Sana et al. (2012) which 
significantly favors short period systems. Such a prefer¬ 
ence had been observed in previo us surveys uncorrected 


for biases (e.g. Mason et al. 2009|), but was generally in- 
terpreted as being the result of selection effects. |Sana| 
et al. (2012) demonstrated that this preference remains 
even after carefully correcting for observational biases. 
We adopt 

/p(logp) oc (logp)’", for logp e [0.15,5.5] (1) 


where p is given in days. For our standard simulation we 
adopt TT = —0.5. We change the slope of orbital period 
distribution from tt = —0.75 in model N-pl to tt = —0.35 
in model N-p2. 

Note that the spectroscopic observations can only re¬ 
liable probe systems with logp 3.5. However, wider 
systems can still produce double compact object merg¬ 
ers as we will show in the following section. For wider 
systems we adopt the simplest assumption we can take 
and extrapolate the distribution, since we have no rea¬ 
sons to believe that the binary fraction suddenly drops 
beyond logp = 3.5. 

This assumption is consistent with the recent find¬ 
ings b y interferometric s tudies of nearby Galactic massive 
stars. Sana et al. (2014) provided a large systematic sur¬ 
vey probing companions of O stars at angular separations 
between 1 and 100 mili-arcseconds. For unevolved mas¬ 
sive stars (O stars with luminosity class V) the detected 
companion fraction reaches 100% at 30 mili-arcseconds. 
The physical separations this corresponds to will remain 
uncertain until more accurate distance measurements be¬ 
come available. Roughly it corresponds to physical sep¬ 
arations of 60-6000 AU. In our preferred units for the 
orbital separation a this corresponds to log a{ Rq) = 4.1- 
6.1. Considering systems with typical masses we find 
that our extrapolation of the orbital period distribution 
and binary fraction are consistent with these observa¬ 
tions. 

Mass ratios — For the distribution of mass ratios, which 
we dehne as the mass of the initially less massive star over 
the mass of the more massive star, i.e. q = M2/Mi, we 
use 

fq{q)(xq'^, forge [0.1,1] (2) 


where k = —0.1 ± 0.6 according Sana et al. (2012). We 
adopt K = 0 such that the distribution becomes uni¬ 
form distribution, which has also been found in several 


recent studies such as 

Kobulnicky & Fryer 

0 

0 

CNl 

) and 

Kobulnicky et al. 

(20U 

t|. To consider the uncertainties 


N-ql to K = —0.7 in model N-q2. 

We note that the most recent observations rule out 
the prese nce of a so-called twin populat ion of equal mass 


systems (Pinsonneault & Stanek 2006|). The idea of the 
possible existence of such a population gained interest as 
it favors the formation of double compact object mergers, 
in particular through the so-called doubl e core formation 
channel proposed by Dewi et al. (2006). However, the 
claimed evidence of such a population has been demon¬ 
strated to be the result of observational biases towards 


equal mass systems that were not accounted fo r (Lucy 
2006 Sana et ^[2013 Cantrell fc Dougan||2014 ). 


Eccentricities — For the eccentricity distri bution in the 
very young open clusters Sana et al. (2012) finds 


/e(e) 


for e G [0.0, 0.9] 


(3) 


with 7] = —0.42 ± 0.17. The very short period systems 
(p < 4 days) show a larger degree of circularization as ex- 
pected from the short time scale for tidal circularization 
( ZahnlloTS ). Unfortunately the data sample is not large 
enough to provide a reliable separation dependent eccen¬ 
tricity distribution. Instead we adopt this distribution 
as initial distribution for our simulations, independently 
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of the period. We explicitly follow the effects of tides in 
our simulations, which quickly circularizes the shortest 
period systems, consistent with the observations. As we 
will show later this assumption is justifiable as the vari¬ 
ations in the eccentricity result in only minimal changes 
in the rates. In our standard simulation N we adopt 
77 = —0.42. We consider uncertainties by changing the 
the slope to rj = —0.59 in model N-el to 77 = —0.25 in 
model N-e2. 


Binary fra ction — The spectroscopic survey by |Sana| 
et al. (20121 yields a binary fraction of /sb = 0.7 ± 0.1 
after correcting for biases. The binary fraction here is 
defined as 

where A^bin and Aijingie are the number of binary systems 
and the number of single stars respectively. 

This fraction refers only to systems that have or¬ 
bital parameters within the considered boundaries , i.e. 
q S [0.1,1], logp(days) S [0.15,3.5] and e G [0.0,0.9]. 
The remainder of the sample 30%) consists of stars of 
unknown nature, including wider binaries, binaries with 
more extreme mass ratios and possibly genuine single 
stars. Mo tivated by the interferometric survey (Sana 
et al.|2014|) we designate th e remainder as wi de binaries. 
We extrapolate the original|Sana et al.|(|2012D period dis¬ 
tribution that extends to logp = 3.5 to logp = 5.5. Such 
extrapolation results in a binary fraction, including wide 
systems, of /bin = 1- We also consider a reduced binary 
fraction of /bin = 0.85 in model N-fl and /bin = 0.7 in 
model N-f2. 

Masses — The o bserved distr i bution of primary masses 
in the sample of Sana et al. (2012) is consistent with 


a standard Kroupa 1 2002 ) mass function. Although we 
only simulate the evolution of massive binaries, for nor¬ 
malization we adopt the three component power law for 
the primary mass or single star ttti, given in solar units 
Mq 


/mi(wi) (X 


1777/^'^, for mi G [0.08,0.5] 
for TTii G [0.5,1.0] 

[ ?n7“, for mi G [1,150] 


(5) 


We adopt a = 2.7 in our s tandard model, consistent with 
the field star p opulation (jKroupa et al.||1993| |Kroupa & 
Weidner 2003). To allow for uncertainties we consider 
a = 3.2 in model N-ml and a = 2.2 in model N-m2. 

2.2. Old standard model (O) and its variations 


(| 2012 ) reference model 
their work systems are 


We adopt the Dominik et al. 
as the old standard model U. fii 
drawn from an initial distribution of separations instead 
of orbital periods. This distribut ion is assumed to be flat 
in the log ( Abt||l98'3 6pik||l9M ) 


This was adopted to ensure that at zero age both stars 
are well within their Roche lobes at closest approach, 
i.e. the periastron distance dper = a(l — e). The distri¬ 
bution of mass ratios is a pseudo-flat distribution. The 
only differe nce with the new dist ribution concerns the 

(|2012|) consider q G [gmin, 1] 
= U.1)8/7711 to avoid selecting secondaries be- 


boundaries. Dominik et al. 


For the distribution of eccentrici- 


with ( 7 m; 

low H-burning limit. 
ties a the rmal-equilibrium distribution is adopted (|Heg-j 
gie||1975 ), which favors eccentric systems: 


/e(e) oc e, for eG [0,1]. ( 8 ) 

The binary fraction adopted by Dominik et al. (2012|) 
is /bin = 0.5, corresponding to the parameter boundaries 
specified above. For comparison we also provide esti¬ 
mates based on the old distribution but assuming a high 
binary fraction of /bin = 1 in model 0 -fl. 

The distribution of primaries masses (and single star 
masses) is identical to the one we adopt in the new ref¬ 
erence model, given in Eq[^ 


3. COMPUTATIONAL METHOD 
For the purpose of this study we use the binary 


lation synthesis code StarTrack (Belczynski et al. 2002 


i 


opu- 


2008), which has been used extensively to simulate the 
evolutionary channels and formation rates of compact 
objects www.syntheticuniverse.org. The suite of our 
simulations is available on this website. 

This code belongs to a family of very fast and stable 
binary evolutionary codes, that provide powerful tools 
to explore the vast multi-dimensional space of the initial 
parameters that determine the fate of a binary system. 
These codes also enable the exploration of the effect of 
model uncertainties. These codes relay on precomputed 
grids of detailed stellar models and approximate treat¬ 
ments of the physical processes as described below. De¬ 
spite the approximations, studies with this and similar 
codes based on the same philosophy have enabled in¬ 
sights into the many exotic phenomena resulting from 
low and high mass interacting binaries (e.g. recent work 


by|Izzard et al. 


niter et al.|2009[|Men- 


nekens et ai.||20101 Vanbeveren et ai. 

|2013| Fragos et al. 

20131 Abate et al. 2013 de Mink et ai. 

2013 20141 1'oonen 

et al. 20141 Kochanek et al. 

20141 Claeys et al.|2014 Men- 

nekens & Vanbeveren|2014 

IRuiter et al.|2014||Schneider 


A full description of the StarTrack code used in this 
study and assumptions concerning the treatment of stel¬ 
lar evolution and binary interact i on ca n be fo u nd in the 


papers abov e Belczynski et al. (2002 2008), Dominik 


et al. (2012) and references therein. Below we provide a 


summary of the main assumptions relevant for this study. 
3.1. Physical assumptions 

This study is a c omparative study of th e impact of the 

is used as our 
012 ), also when 


initial conditions. Dominik et al 


reference model. We use Dominik 


/a(loga) oc constant, for logo G [aniin)5]. (6) assumptions are uncertain. 


lioTl 

kal.|(| 2 i 


The lower boundary Oniin is assumed to be a function of 
the stellar radii of the primary and secondary star at the 
zero age main-sequence (ZAMS), 


Stellar evolution — The code is based on detailed, but 
non rotating , single stellar evolutionary models by |Pols| 


et al. (1998). For this purpose the code utilizes the tit 


Q-min — 


_ 2i?zAMS.l + 2i?zAMS,2 


formula to the detailed models by Hurley et al. (|2000| 


1 — e 


(7) with several adaptations described m Belczynski et al. 
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(2002). The effects of mass loss and accretion are sim¬ 
ulated using algorithms origi nally developed by (Tout, 
et al.|19^ Hurley et al.|2002 ) to account for effects such 
as rejuvenation ot the accreting star when appropriate. 

Stellar wind mass loss — We employ updated wind mass 
loss ra tes, which include winds mass loss fro m early-type 
stars (Vink et a l. 2001), Wolf-Rayet sta rs (Hamann & 
Koesterke 1998 ' Vink & de Koter 2005) and enhanced 
mass loss rates for Luminous Hlue Variables calibrated 
to account for the observed distribution of masses of 


265 km s based on the observed distributions of ra¬ 


know n black holes (Belczynski et al. |2010 Orosz et al. 
2011). As a result the simulations allow tor the torma- 
tion of black holes with masses up to 15 Mq in a solar 
metallicity environment. For sub-solar metallicity en¬ 
vironments {Z = 0.006), black holes up to 30 Mq are 
produced, consistent with the mass of the most massive 
known stellar-origin black hole in the ICIO X-1 system 
(Prestwich et al. 2007 Silverman & Filippenko |2008). 
Under these assumptions the simulations predict the tor- 
mation of black holes with masses up to 80 Mq in low 
metallicity {Z = 0.0002) environments for stars below 
initial mass of 15OM0. 


Roche-lobe overflow and common envelope evolution — To 
determine whether mass transfer through Roche-lobe 
overflow is stable we consider t he stellar type of the do nor 
and mass ratio as outlined in (Belczynski et al.|2008 ). In 
the case of stable mass transfer, we assume that half 
of the mass lost by the donor is accreted by the com¬ 
panion. The remainder is leaving the system carrying 
a specific angular momentum of 27ra^/P, where a and 
P are the orbita l separ ation and period following |Pod- 
siadlowski et al. (1992). Common envelope evolution is 
acco unted tor using the classical energy balance formal¬ 
ism (|Webbink||1984D adopting a value of acE = 1 for 
the envelope efficiency parameter. The parameter A de- 
scribin g the binding energ y of the envelope is taken from 
fits by jXu fc Lij ( |2010a 0, implemented as described in 


Dominik et al. '( (2012[ ). 


Remnant masses for compact objects — To obtain the final 
mass of the compact objects, we consider the formation of 
a proto-neutron star and subsequent fall back of material. 
These are inferred to be a function of the carbon oxygen 
core mass as outlined in Section 4.2 of Fryer et al. (2012). 
The employed scheme utilizes the rapid supernova model 
that is able to reproduce the apparent mass gap between 
NSs and BHs ( Belczynski et al.|2012b ). We also account 
for the formati on ot neutron stars by electron capture 
supernovae (e.g. Nomoto|1984 iPodsiadlows k i et al . 2004) 
4.4 oiil''rver et al. (12012). 

fo 


as outlined in Section 4.4 oi|Fryer et aT 


For 

single stars this effectively results in the formation of 
neutron stars for initial stellar masses above 7.6 Mq. The 
transition from neutron star to black hole is set at a 
mass of 2.5 Mq for the resulting compact object. For 
high metallicity {Z = 0.02) this makes single stars above 
initial mass of 2IM0 to form BHs. 

Supernova ki cks — We account for the clas sical Boersma- 
Blaauw kick (Boersma 1961| Blaauw 1961) resulting from 
mass loss during the explosion. In addition, newly 
born compact objects receive natal kicks which are given 
by a Maxwellian distribution with a ID rms of cr = 


dio pulsars (Hobbs et al. 2005). The kicks are lowered 
proportional to the amount ot f all back as described in 
Section 4.5 of Fryer et al. (2012). An exception is made 
for neutron stars formed through electron capture super¬ 
novae for which no natal kicks are assumed in this simu¬ 
lation. As a result of this the mergers of double neutron 
stars are biased toward neutron stars formed through 
electron capture supernovae. 

Submodel A and B — To consider the large uncertainties 
concerning binaries that begin Roche-lobe overflow while 
the donor star crosses the Hertzsprung gap, we consider 
two extreme cases. In submodel A we consider all com¬ 
mon envelope formation channels including donor stars in 
any post main sequence evolutionary stage. In submodel 
B we explicitly exclude all evolutionary channels where 
the common envelope is initiated by a Hertzsprung gap 
donor star. Submodel B can be considered as a more con¬ 
servative estimate of t he merger rates for NS- NS, BH-NS, 
and BH-BH binaries ( Belczynski et al.|2007 ). A detailed 
study and revision of criteria tor the onset 01 common en¬ 
velope will be presented shortly (Belczynski, Pavlovskii 
& Ivanova, in prep.). 

Metallieity — With the reach of advanced LIGO/Virgo 
mergers will originate from all sorts of environments with 
high and low metallicity. We therefore consider two 
representative metallicities Z = 0.02 (typical of solar 
neighborhood) and Z = 0.002 (typical of small starburst 
galaxies). The physically consistent way to derive merger 
rates would be to start from an appropriate star forma¬ 
tion history in the Universe, adopt a model of metallicity 
evolution with redshift and then calculate broad range of 
population synthesis models with varyi ng metallicity in 
order to obtain current rate of mergers (jO’S haughnessy 
et al.|[2008 2010 Dominik et al.|[2013 


2015). Given the 


scope ot this study we adopt a sirnpler approach and take 
an even mix of these two compositions as a crude rep- 
rese ntation of metallici ty distribution in local Universe 
(e.g. Banter et al.|2008 ). Although simplified, this is suf- 
ficient to judge the overall impact of the change in initial 
conditions and their associated uncertainties. 


Difference with respect to Dominik et al. (2012) — Since 
the publication of t he po pulation synthesis calculations 
by Dominik et al. (2012), it was found that a techni¬ 
cal Dug~waFln7roauced~during one of the annual code 
updates, related to the treatment of tidal locking. For 
detached binaries, in which both stars reside well within 
their Roche lobe, the effect of tides on the stellar and 
orbital spin is negligible. In such systems, the spins of 
the stars are typically not synchronized with the orbital 
motion. As the stars evolve they change their spin pe¬ 
riod, for example as a result of evolutionary expansion. 
At some point the spin and orbital period may become 
comparable. Such a system appears to be synchronized, 
at least momentarily, eve n though tides are no t effec¬ 
tive. The simulations by Dominik et al. (2012) incor¬ 
rectly treated these systems as tidally locked, resulting 
in an incorrect further orbital evolution for specific cases. 
In general, this resulted in shrinking of the binary sys¬ 
tem with the overall effect of over predicting the double 
compact object merger rates. 
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For majority of the cases the differences are small: 
within factor of ~ 2, i.e. comparable to differences aris¬ 
ing from the change of initial conditions. This can be 
seen when comparing the r esults before th e corre ction 
given in Table 2 and 3 of Dominik et al. (2012), the 
standard model marked ”S’ m their paper and our re¬ 
sults for the same assumptions, but after correction for 
the bug, as we provide in our Table 1, marked as the old 
standard model, ”0”. 

In two cases, both for solar metallicity submodels 
B, involving binary merges with black holes the dif¬ 
ferences are larger. The BH- NS Galactic merger r ate 
was revised from 0.2 Myr“^(Dominik et al. 2012) to 
0.06 Myr“^(current study). The BH-BH Galactic 


merg er rate was revised from 1.9 Myr (Dominik et al. 
2012) to 0.22 Myr“^(current study). 

Although the changes for these high metallicity chan¬ 
nels are substantial, we repeat that the overall rate is 
completely dominated by the low metallicity channels. 
Therefore the impact of this correction on the overall 
cor npact merger rates, masses and de lay times presented 
by ( Dominik et al.|[2012 2013 2015) is negligible. 


3.2. Computation of the merger rates 

We randomly draw binary systems from the initial dis¬ 
tribution functions, once for high metallicity {Z = 0.02) 
and once for low metallicity [Z = 0.002). For com¬ 
putational efficiency we only select binary systems that 
are massive enough, within some safety margin, to po¬ 
tentially produce a double compact object involving a 
neutron star or a black hole. For primaries we take 
Ml > 5 Mq and for secondaries M 2 > 3 Mq. 

Using the StarTrack population synthesis code we 
evolve Nsira = 2 X 10® of such binary systems. We esti¬ 
mate the mass formed in our simulations. To obtain the 
total mass in stars down to the lower stellar mass limits 
we integrate over the full extent of the initial mass func¬ 
tion (0.08— I5OM0) and mass ratio distribution. We in¬ 
clude the appropriate mass for single stars when needed. 
For this we assume that for each binary system with a 
primary mass Mi, there are (1 — /bin)//bin single stars 
with a mass M = Mi 

We present our rates in terms of the rate for a fiducial 
G alaxy, . We chose t o adopt the exact same method 
as Dominik et al. (2012) to allow for comparison. We 
therefore assume a lOGyr continuous star formation rate 
^SFR = 3.5 Mq Myr“^. The resulting total mass formed 
in stars is within a factor of two of the estimates for 
the present-day stellar mass in the Milky WajH The 
star formation h istory of our galaxy is not wellknown 

(e-g- 

epoc 


Wyse 2009) but there is evidence for s everal discrete 
IS of enhanced of star formation (e.g. Gignoni et al. 


20061. The assumption of continuous star formation may 
be a rather crude approximation for the Milky Way itself. 
It may however serve as the rate for a fiducial Milky 
Way-like galaxy obtained after averaging over the LIGO 
detection volume. 

In Table we provide the rates for the different simu¬ 
lations. These rates can be converted into approximate 
volumetric rates, as used by the LIGO/Virgo collabora- 


^ Estimates for the total stellar mass by|Mjmnet^am2006| yields 
4.85-5.5x IQio Mq and 6.43±0.63 X lO^o Mq by |McMillan| 12011 1. 


tion, with the following expression. 


'R-woX = 10 yr ^ Gpc ® 


p 

TT-gal 

.O.OlMpc”®. 

-1 

7 


(9) 


where p is the local density of Milky Way like galaxies 
and i?gai the Galactic merger rate. 

The detailed choices on how to normalize the rates 
are, to some extent, arbitrary. This is a result of the 
uncertainties in the low mass binary fraction and the 
poorly constrained distribution of their mass ratios. We 
chose to adopt very simple assumptions and we assume 
for the sake of this study that both binary fraction and 
mass ratio distribution adopted in a given model for mas¬ 
sive stars hold true for the entire considered mass range 
(0.08 - 150 Mq). 

Our Monte Garlo simulations are subject to statisti¬ 
cal fluctuations due to they finite size. We estimate this 
by counting the total number of double compact objects 
that merge with 10 Gyr in each simulation, iVx, consider¬ 
ing NS-NS, BH-NS, BH-BH separately. We use l/y/Nf 
as an approximation for the statistical uncertainty. 


4. RESULTS: DISTRIBUTIONS 
4.1. Comparison of the initial distributions 

The effective initial distributions for all systems with 
a primary mass of at least 7 Mq are shown in Fig. for 
the old and new assumptions. Systems with lower pri¬ 
mary masses are not capable of producing double neutron 
stars in our simulations. The histograms are shown for a 
metallicity Z = 0.02 and normalized to unity. The main 
differences between the old and new assumption concern 
(a) the period/separation distribution and (b) the eccen¬ 
tricity distribution. 

Periods and separations — The main difference between 
old and new assumptions concerns the tightest systems. 
The new distributions show a clear preference for short 
periods / small separations as can be seen in central top 
and bottom panels of Fig. The old distribution is flat 
for large separations but tight systems (a < 100 i?© ) 
are strongly suppressed. This is a result of the adopted 
lower limit Omin, which depends on the eccentricity and 
the stellar radii ( see Sect. [^. Such a turnover is not 
observed for massive binaries. Even though we do not 
know at present how such close systems form, they exist 
and must be accounted for in the simulations as we do 
in the new initial conditions. 

Eccentricities — The second main difference concerns 
the eccentricities (depicted in the top right panel of 
Fig. 04 The distribution adopted in the old simulations 
is the thermal distribution which strongly favors eccen¬ 
tric systems. The new distribution favors (near) circular 
systems. The high eccentricities in the old simulations 
allow much wider systems to still interact as discussed 
below. 

Periastron separations — The distribution of periastron 
separations dper = ct{l — e), i.e. the distance of closest 
approach at zero age, is shown in the bottom right panel 
in Fig.[2 To first order this is the most relevant parame¬ 
ter governing when systems start to interact by tides and 
later mass transfer. The old simulations favor wider and 
more eccentric systems, while the new simulations favor 
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Fig. 1. — Distributions of initial parameters for the simulated systems resulting from the old and new assumptions. All histograms are 
normalized to unity. 


tighter and more circular systems. By coincidence, these 
changes partially compensate each other. 

The old assumptions resulted in a nearly flat dis¬ 
tribution of periastron separations extending from 
logdper(R©) = 1-5, being nearly flat between 

logdper(R.0) = 1.5-4 and turning over on both ends. 
Instead the new distributions peak at short periastron 
separations of about 1O-5OR0 and gradually decay for 
higher separations. 

However, in the regime of relevance for the formation 
of double compact objects (logdper(R0) = 1.5-4, as we 
will see in the next section) both, old and new, peri¬ 
astron distribution are rather similar. We note only a 
slight domination of old distribution in this regime. The 
distributions in this figure are normalized to unity. The 
larger overall binary fraction in the new simulations will 
compensate for this. 

Primary masses and mass ratios — In both simulations 
the mass ratios are drawn from a flat seed distribution 
(the bottom left panel of Fig. f^. They only differ in the 
adopted boundaries. In the md simulation companion 
masses are drawn down to the hydrogen burning limit. 
In the new simulations the minimum mass ratio adopted 
is 9min = 0.1. In our simulations with old initial distri¬ 
butions we found no double compact objects resulting 
from systems with such extreme mass ratio. This differ¬ 
ence therefore only has a very minor effect on the final 
normalization. 

The distributions of initial primary masses Mi (shown 
in the top left panel of Fig. [I) are identical in both sim¬ 
ulations. The only differences are stochastical in nature. 

4.2. Comparison of the birth properties of the 
progenitors of double compact object mergers 

Only a very small subset of the simulated binaries pro¬ 
duces two compact objects that are bound and that will 
merge within 10 C^r. The birth properties of this subset 
are shown in Fig. where we compare the impact of the 


old and new assumptions for the initial conditions. 

Note that the histograms are normalized to unity to 
allow comparison of the shape of the distribution, the 
resulting rates are discussed below. We separately show 
results for low metallicity (Z = 0.002, top panels) and 
high metallicity (Z = 0.02, bottom panels). 

Since the trends observed when comparing the old and 
new initial distributions are very similar for submodel 
A and B, we chose to only show the results submodel 
B here. This is the more conservative submodel which 
excludes progenitor systems in which a Hertzsprung gap 
star acts as donor star during CE phase, see Section 3.1. 
In practice this means that on average the progenitors in 
model B result from wider systems and show a slightly 
stronger preference for equal mass progenitor systems. 
Data for submodel A (as well as data subdivided by 
the double compact object type) can be obtained from 
the online data repository for readers that are interested. 
More insight in the differences between submodel A and 
B can be found in (Dominik et al.|[2012|. 


Primary masses — The distribution of initial primary 
masses of the progenitors consists of separate compo¬ 
nents. The narrow low mass peak around IOMq is gener¬ 
ated by the progenitors of double neutron stars. The high 
mass component is only prominent in the low metallicity 
simulations. This component consists of the progenitors 
of mergers that involve a black hole. 

The minimum initial mass to form a black hole is 4OM0 
for the primary star of a double compact object progeni¬ 
tor at Z — 0.02. This is much higher than the minimum 
mass for a single star (or star in a very wide non inter¬ 
acting binary) to form a black hole in our simulations, 
which is 2IM0 at Z = 0.02. The difference is the result 
of the additional mass loss experienced in binary sys¬ 
tems due to Roche lobe overflow and common envelope 
ejection from BH progenitor. At lower metallicity these 
minimum masses decrease. At Z = 0.002 we find a min¬ 
imum initial mass of about 2OM0 for the primary star 
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Fig. 2.— Birth distributions of initial binary parameters for the progenitors of double compact object mergers. The normalized histograms 
show the distributions for submodel B combining all types of mergers, for low and high metallicity. 


of a double compact object progenitor involving a black 
hole. 

The difference between the old and new simulations 
is very small. At low metallicity there is a slight ten¬ 
dency towards lower mass progenitors (favoring NS-NS 
progenitors) in the new simulations. 

Mass ratios — For both the old and new simulations we 
find that the progenitors exclusively result from systems 
with mass ratios larger than about M2/M1 > 0.4 at low 
metallicity and about 0.6 at high metallicity. For both 
metallicities there is a preference for progenitors with 
very equal masses, i.e. 0.9 and higher. 


The difference between the old and new simulations is 
marginal, unsurprisingly given the similarity of the input 
distribution. A very slight tendency for the old simula¬ 
tions to favor more equal mass systems is observed at 
high metallicity, while the opposite behavior is observed 
at low metallicity. However, the differences are insignih- 
cant in comparison with the model uncertainties. 

Periods, separations and eccentricities — The initial pe¬ 
riod and separation distributions of the progenitor sys¬ 
tems show a clear difference between the old and new 
initial conditions. For both metallicities the old simu¬ 
lations are biased towards wider systems than the new 
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Fig. 3.— Variations in the merger rate predictions double neutron stars for the old and new standard assumptions and the effect of the 
allowed variations in the initial distributions of binary parameters. See also Table 2. Note that the Poisson variations (indicated as error 
bars) resulting from the finite size the number of simulated systems are comparable with some of the variations. We use a double mirrored 
logarithmic vertical axis to visualize the very small variations we obtain in most cases. For comparison, uncertainties from other sources 
(e.g., evolutionary uncertainties) are typically an order of magnitude or more. 


simulations. At high metallicity the distributions peak 
near 15OR0 in the new simulations and near 300Rq in 
the old simulations. At low metallicity the distribution is 
clearly bimodal peaking around 800 and SOOOR© in the 
new simulations. 

The distribution of initial eccentricities of merger pro¬ 
genitors closely reflects the shape of the input distribu¬ 
tion of the population, although in both cases systems 
with larger eccentricities are over represented. 

Periastron separations — Although the distributions of 
initial separations a and eccentricities e are very different 
between the old and new simulations, the distribution of 
initial periastron separations dper of merger progenitors 
is nearly indistinguishable. 

For low metallicity the distribution is again bimodal. 
Detailed inspection of the evolutionary paths of the sys¬ 
tems that populate both peaks show that the tighter sys¬ 
tems first evolve through a stable mass transfer phase. 
The wider systems first evolve through a common enve¬ 
lope phase. These evolution ary channels are ident ified 
and described in earlier work (Dominik et al. (20121; see 
their Table 5 and associated text). 

Both channels exist in the old and new simulations, 
since the input physics assumptions are not changed. 
The similarity in the periastron separation distributions 
implies that the relative contribution of both channels 
did not significantly change when switching from the old 
to new assumptions. 

4.3. Distributions of final properties 


The distribution of the final properties predicted by 
our simulations that are of interest for gravitational wave 
searches include the delay times, component masses, 
mass ratios and chirp masses. Detailed models of the 
anticipated waveforms will enable the infere nce on sev¬ 
eral parameters, such as component masses (Aasi et al. 
2013a| IVeitch et al.||2015|) and if multiple events are de¬ 
tected const raints on th^roperties of the popu lation can 
be inferred ( Mandel||20ici[ Mandel et al.||2015 ). 

We find tnat the distributions of hnaf properties are 
remarkably insensitive to the assumed initial conditions. 
The predicted distribution of the component masses, 
mass ratios, chirp masses and delay times show no sig¬ 
nificant changes when the initial conditions are varied 
within current uncertainties. They are alm ost identical 


to those presented by Dominik et al. (2012), to which we 


refer for detailed discussion. 

The insensitivity of rate predictions results from the 
(coincidental) similarity of the initial periastron separa¬ 
tion distribution and the mass ratio distribution, which 
remained practically unchanged. These two parameters 
are key in determining the fate of the binary system (and 
thus merger rates). However, the physical properties of 
merging binaries are not sensitive to changes in initial 
conditions. 

There are relatively few formation channels for double 
compact object mergers. For a given model of evolution, 
these formation channels originate from a rather narrow 
volume of initial parameter space (“the formation vol¬ 
ume” ). Change of initial conditions alters the number of 
the progenitor binaries in the formation volume, but has 
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Fig. 4.— As Fig. I^for the merger rates BH-NS (left) and BH-BH (right) systems. 


very little effect on the binary merger properties. This 
insensitivity is not just connected to the specific initial 
conditions considered here, but is of more fundamental 
nature. It is a robust result that has also been found 
in connec tion to other compact object binaries. For ex¬ 
ample by Kalogera & Webbink (1998) in the context of 
low-mass X-ray binaries, and see also |Belczynski et al.| 
([ 200 ^. 

The robustness of the final distributions against initial 
condition uncertainties considered here, but also wider 
variations, is reassuring in the light of anticipated grav¬ 
itational wave detections. It removes one layer of un¬ 
certainties and allows for more optimism concerning the 
potential of such detections for constraining more inter¬ 
esting aspects such as the physics uncertainties. 

The final properties of our simulations are avail¬ 
able through our online data files available at 
(www.syntheticuniverse.org). 


5. RESULTS: MERGER RATES 

In Table we list the double compact object merger 
rates for our new standard model N, the old standard 
model O, and variations on the models that explore the 
effect of uncertainties in initial distributions. All rates 
are expressed as th e ra te for a fiducial Milky Way Galaxy 
as detailed in Sect |3.2[ We give the results both for sub¬ 
model A and the more conservative submodel B. Merger 
rates are also presented separately for low {Z = 0.002) 
and high metallicity {Z = 0.02) stellar and binary evo¬ 
lution. 

In Table we present merger rates that correspond to 
a fiducial Galaxy with an even mix (50%-50%) of low- 
metallicity {Z = 0.002) and high-metallicity (Z = 0.02) 
stars. Absolute mixed merger rates are given only for 
our new standard model (an average of low- and high- 
metallicity fiducial Galactic rates from Table and can 
be readily obtained for any other model. For all other 
models only relative change of rates with respect to our 
new standard is given. A graphical representation of the 
changes relative to the new standard model is given in 
Figure 


The overall conclusions are: (i) merger rates increase 
slightly with the new assumptions, and they are within a 
factor of 2 of old standard rates; (ii) merger rates change 
by less than factor of 2 within associated uncertainties 
embedded in new initial distributions concerning binary 
parameters (the only exception being the IMF). These 
conclusions hold for all types of double compact object 
mergers (NS-NS, BH-NS and BH-BH) and the changes 
are typically much smaller than the maxima listed above. 
In the following section we discuss the various trends. 


5.1. Impact of the new initial conditions 

The adoption of the new assumptions for the initial 
conditions leads to higher estimates for all merger rates. 
The largest change is found when adopting model A. The 
old rates are 39% lower than the new rates for NS-NS 
mergers, 29% lower for BH-NS mergers and 27% lower for 
BH-BH mergers. An overview of the relative changes is 
given in Tab. [fusing the new standard rate as reference. 

When adopting submodel B the observed changes are 
even smaller in the case of BH-NS and BH-BH merg¬ 
ers. The submodel B excludes progenitor channels which 
involve com mon envelope evolution with a Hertzsprung 
gap donor (Belczynski et al. 2007). These typically 
come from progenitor binaries with smaller initial sepa¬ 
rations where the old and new distributions deviate most 
strongly. Excluding these progenitors reduces the differ¬ 
ence between old and new simulations. Eor the BH-NS 
in submodel B the changes are so small that they are 
of the same order as the expected statistical fluctuations 
resulting from the finite number of systems in our simu¬ 
lations. These average statistical fluctuations are shown 
in Eigure|^as an errorbars. 

The main reason for the slight increase in the merger 
rate is the large binary fraction adopted in the new ini¬ 
tial conditions {ft = 1). The old distributions with the 
high binary fraction of /(, = 1 (model 0-fl) give merger 
rates that are higher (up to ~ 50% in the case of BH-NS 
mergers in submodel B) than for the new standard model 
(see Table 2). 
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5.2. Variations due to uncertainties in the new initial 
conditions 

All relative merger rate changes are smaller than 34% 
for all types of double compact objects when varying 
the initial conditions concerning binary properties (mod¬ 
els Nml, N-p2, N-ql, N-q2, N-el, N-e2, N-fl, N-f2; see 
Tab.[^. This is not true for variations in the IMF, which 
we wul discuss at the end of this paragraph. 

The initial distribution of mass ratios that rather 
steeply falls off with mass ratio (model N-q2) is the cause 
of merger rate decrease of the order of ~ 30% for NS-NS 
and BH-BH mergers. A distribution favoring equal mass 
systems at birth mildly favors the formation of NS-NS 
systems and BH-BH systems. The BH-NS mergers show 
the same trend, but the changes are comparable to the 
Poison fluctuations. Varying the initial period distribu¬ 
tion has a smaller effect, below 20 % and approaching 
our statistical uncertainties. A flatter period distribu¬ 
tions, i.e. one that favors wider systems (model N-p2), 
results in a more effective formation of compact object 
mergers. A period distribution favoring short period sys¬ 
tems more strongly (model N-pl) mainly adds systems 
that interact or merge prematurely, preventing the for¬ 
mation of a double compact object system. Varyin g the 
eccen t ricity distribution within the ranges given by [Sana 
et al. ( 2012 ) has no significant effect on the merger rates. 
All variations found are less than 5% for all submodels 
and all types of mergers considered. These variations are 
dominated by the statistical fluctuations in our Monte 
Carlo simulations. A reduction of the total binary frac¬ 
tion to /bin = 0.85 (which roughly corresponds to the 
one sigma lower limit of intrinsic sp ectroscopic binary 
fraction, /gp = 0.7 ± 0.1, derived by |Sana et al.||2012| 
leads to a reduction of all compact object merger rates 
by 10% (model N-fl). A further reduction of the total 
binary fraction to /bin = 0.7 (model N-f2), which roughly 
corresponds two the two sigma lower limit on the obser¬ 
vational constraints, reduces the merger rates by 15%. 

Significant change of the power law exponent of the 
IMF for stars more massive than IMq generates the 
largest change of double compact merger rates in our 
suite of models. We allow for a change of the exponent 
by 0.5 up and down from our standard value a = 2.7. 
Note that this affects the mass distribution of the pri¬ 
mary stars in binary systems and single stars for those 
models where /bin < 1 (model 0-fl, N-fl, N-f2). Merger 
rates decrease within a factor of 2.7 for NS-NS mergers, 
a factor of 4.8 for BH-NS mergers and a factor of 5.8 
for BH-BH mergers for our most steep IMF (model N- 
ml) in respect to our new standard model (N). Merger 
rates increase within a factor of 2.2 for NS-NS mergers, 
a factor of 4.1 for BH-NS mergers and a factor of 4.7 for 
BH-BH mergers for our most flat IMF (model N-m2). 
These changes are due to two effects. First, the steeper 
the IMF the fewer double compact objects form in a pop¬ 
ulation of binary stars. The heaviest BH-BH mergers are 
affected the most. Second, the slope of high-mass IMF 
has a significant impact on the total mass of the fiducial 
Galaxy. For the same number of massive binary stars 
that we draw from initial distributions {Nsim = 2 x 10 ®) 
the corresponding total mass of simulation (over the en¬ 
tire stellar mass range) is ^ 2.5 higher and ~ 2.5 smaller 
for our model with the steep and flat IMF, respectively 


as compared with our adopted standard IMF. 

6 . DISCUSSION AND CONCLUSION 


Recently, Sana et al. (2012) have presented measure¬ 
ments of iiiitial binary parameters for massive stars. 
These measurements, in some aspects, are in stark con¬ 


ies of massive stars (e.g. Belczynski et al. 2002 Voss & 

Tauris 20031 

Dominik et al. 2012 Mennekens & Van- 

beveren 2014 

). The old standard prescriptions (a hat 


tribution strongly favoring eccentric systems and typi¬ 
cally an adopted binary fraction of 50% counting sys¬ 
tems with separations up to 10® Rq) may still apply for 
low- and intermediate-mass stars, but the new measure¬ 
ments show very different initial binary parameter distri¬ 
butions for massive stars (strong overabundance of very 
tight systems, a binary fraction that reaches 100 % when 
including systems up to about lO'^'^R© and mostly cir¬ 
cular systems). 

We find a small increase in the NS-NS, BH-NS and 
BH-BH merger rate as a result of the new initial condi¬ 
tions (due to the increased binary fraction), with the old 
rates being up to 40% lower than the new rates. Sec¬ 
ondly, we study the impact of the allowed variations of 
the new binary parameter distributions within the obser¬ 
vational uncertainties. We find that they have negligible 
impact (typical variations of a few percent, always less 
than a factor of 2). All changes are negligible in compari¬ 
son with the evolutionary model uncertai nties (which are 
typic ally 1-2 orders of magnitude; e.g., [Dominik et al.| 


2012 ). 

The only exception to our above conclusion are the 
variations of double compact merger rates due to the 
uncertainty in the slope of the IMF for stars more massive 
than 1 M 0 . These variations can be as high as factor 
of ^ 3,4 ,6 up and down for NS-NS, BH-NS and BH- 
BH mergers, respectively from our new reference model 
that employs the best estimates on the initial conditions 
for massive binaries. This is due to the change of the 
overall mass of the stellar mass content which is used 
to normalize the rates and to the diminishing number of 
double compact object progenitors with the increasing 
steepness of the IMF. 

The small impact of the new distributions of binary 
properties on the merger rates is somewhat counter¬ 
intuitive. This result stems from the fact that major 
changes in the initial binary parameter distributions can¬ 
cel out in the regime important for close double compact 
object formation. Although, new distributions show a 
much stronger bias toward close binary orbits in compar¬ 
ison to the old distributions, the periastron distance dis¬ 
tribution is similar to the one generated by old distribu¬ 
tions in the regime 1.5 < log(dper/R 0 ) < 4 (see Fig. [^. 
This is because of the strong preference for very eccen¬ 
tric systems adopted in the older distributions, for which 
there is no supported in the new observations (mostly 
circular orbits). This regime of periastron separations 
plays a crucial role in the formation of stellar-mass dou¬ 
ble compact object mergers (see Fig. 1^. The periastron 
distance is key in determining when tides become impor¬ 
tant and when the primary star fills its Roche lobe for the 
first time. Since the distribution of component masses is 
very similar, it is the periastron separation that largely 
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TABLE 1 

Merger rates eor a eiduoial Milky Way Galaxy [Myr“^] 


Z = 0.002 Z = 0.02 


Model ID 

NS-NS 

BH-NS 

BH-BH 

NS-NS 

BH-NS 

BH-BH 

Description 


A 

B 

A 

B 

A 

B 

A 

B 

A 

B 

A 

B 

N 

12.8 

3.68 

4.45 

2.32 

98.3 

15.6 

36.9 

12.5 

1.77 

0.06 

9.73 

0.22 

new standard 

O 

8.17 

2.43 

3.17 

2.24 

72.2 

13.3 

22.3 

7.82 

1.24 

0.01 

6.98 

0.22 

old standard 

O-fl 

13.2 

3.92 

5.11 

3.60 

116 

21.5 

35.9 

12.6 

1.99 

0.02 

11.3 

0.36 

old, /bin = 1-0 

N-pl 

12.5 

3.34 

4.72 

2.23 

90.8 

13.3 

35.9 

11.6 

1.46 

0.05 

9.46 

0.29 

new, (logp)”^’^^ 

N-p2 

13.4 

3.67 

4.64 

2.75 

102 

17.5 

37.3 

12.6 

2.00 

0.05 

9.87 

0.20 

new, 

N-ql 

16.3 

4.28 

5.15 

2.53 

113 

18.8 

44.7 

15.6 

1.48 

0.07 

11.7 

0.24 

new, 

N-q2 

8.65 

2.52 

4.00 

2.23 

69.8 

11.1 

25.7 

8.17 

1.32 

0.06 

6.84 

0.13 

new, q ' 

N-el 

13.2 

3.71 

4.49 

2.40 

99.0 

16.0 

38.8 

12.9 

1.70 

0.06 

9.93 

0.22 

new, e 

N-e2 

12.7 

3.41 

4.66 

2.37 

95.3 

15.7 

35.5 

12.0 

1.57 

0.04 

9.90 

0.12 

new, 

N-fl 

11.6 

3.34 

4.03 

2.09 

89.1 

14.1 

33.4 

11.4 

1.60 

0.06 

8.84 

0.20 

new, /bin = 0.85 

N-f2 

10.9 

3.13 

3.78 

1.97 

83.4 

13.2 

31.3 

10.6 

1.51 

0.05 

8.26 

0.19 

new, /bin = 0.7 

N-ml 

4.70 

1.43 

0.99 

0.49 

17.2 

2.77 

13.4 

4.65 

0.31 

0.01 

1.54 

0.04 

new, 

N-m2 

26.5 

7.80 

17.9 

9.51 

457 

66.9 

81.3 

26.9 

7.86 

0.32 

51.4 

1.19 

new, 


Note. — ^Rates correspond to a galaxy with 10 Gyr of constant star formation rate at the level of 
3.5M0 yr“^ as detailed in Section 3.2. Rates are given for low and high metallicity and for two different 
model assumptions. In model A we allow for high rate of common envelope survival, while more the 
conservative model B does not allow the formation of double compact objects for donors on Hertzsprung 
gap. 


TABLE 2 

Relative Changes oe Mixed-Metallicity Merger Rates'^ 


Model ID 

NS 

NS 

BH 

-NS 

BH 

BH 

Description 


A 

B 

A 

B 

A 



N 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

new standard 


(24.85) 

(8.090) 

(3.110) 

(1.190) 

(54.02) 

(7.910) 

absolute rates [Myr“^] 

O 

0.613 

0.633 

0.709 

0.945 

0.733 

0.855 

old standard 

O-fl 

0.988 

1.021 

1.141 

1.521 

1.178 

1.382 

old, /bin = 1-0 

N-pl 

0.974 

0.923 

0.994 

0.958 

0.928 

0.859 

new, (logp)”®'^^ 

N-p2 

1.020 

1.006 

1.068 

1.176 

1.036 

1.119 

new, (fogp)"*^-®® 

N-ql 

1.227 

1.229 

1.066 

1.092 

1.154 

1.204 

new, q^'^ 

N-q2 

0.691 

0.661 

0.855 

0.962 

0.709 

0.710 

new, q ' 

N-el 

1.046 

1.027 

0.995 

1.034 

1.008 

1.025 

new, 

N-e2 

0.970 

0.952 

1.002 

1.013 

0.974 

1.000 

new, 

N-fl 

0.905 

0.911 

0.905 

0.903 

0.907 

0.904 

new, /bin = 0.85 

N-f2 

0.849 

0.849 

0.850 

0.849 

0.848 

0.846 

new, /bin = 0.7 

N-ml 

0.364 

0.376 

0.209 

0.210 

0.173 

0.178 

new, rrLi~^‘^ 

N-m2 

2.169 

2.145 

4.141 

4.130 

4.706 

4.304 

new, rrLi~^‘^ 


Note. — Rates correspond to a fiducial Milky Way Galaxy with an even mix (50%-50%) of 
low-metallicity {Z = 0.002) and high-metalliciW {Z — 0.02) stars (the average of low- and high- 
metallicity fiducial Galactic rates from Table [R. We give all rates relative to our new standard 
model. A graphical representation is shown in~Figure[^ 
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determines the future evolution of the binary system. It 
happens (by pure coincidence) that the old prescriptions 
and the new measurements result in rather similar or¬ 
bital separations at the time of the first mass transfer 
for progenitors of double compact objects. Additionally, 
the observational uncertainties of the new estimates of 
initial binary parameters for massive stars are not large 
enough to cause significant NS-NS, BH-NS and BH-BH 
merger rate changes. Similarly, we find that the distribu¬ 
tions of final properties (delay-times and final component 
masses) are insensitive to the variations of initial condi¬ 
tions. The summary of our merger rate analysis is given 
in Fig. m 

Our study allows to eliminate one layer of uncertainties 
involved in population synthesis predictions for double 
compact object mergers. Such predictions are infamous 
for their number of poorly constrained model parameters. 
Recent years of intensive studies of double compact ob¬ 
ject formation have proven that it is extremely hard (yet, 
not impossible) to identify a well posed non-degenerate 
problem with population synthesis. This fact does not 
reflect the intrinsic weakness of the population synthesis 
method, but rather the lack of understanding of some ba¬ 
sic processes involved in the evolution of massive single 
(e.g., convection, rotation, supernova/core collapse) and 
binary stars (e.g., tidal interactions, common envelope 
evolution). 

The reduction of the modeling uncertainties, like the 
one presented here, is crucial for advancement of our un¬ 
derstanding of the remaining poorly constrained physi¬ 
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